'''
# 10/3 GPT-40-Mini 

import numpy as np
from scipy.integrate import solve_ivp
from lmfit import Model, Parameters
import matplotlib.pyplot as plt

def pk_pd_model(t, y, Ka, Ke, V, Tlag, a, b, E0, g, AaDose):
    Aa, A1, E = y
    dAdadt = -Ka * Aa
    dA1dt = Ka * Aa - Ke * A1
    if t >= Tlag:
        dEdt = E * (a - b * np.log(E / E0)) - g * (A1 / V) * E
    else:
        dEdt = E * (a - b * np.log(E / E0))
    return [dAdadt, dA1dt, dEdt]

# Define the model parameters (initial guesses)
params = Parameters()
params.add('Ka', value=4.0488)
params.add('Ke', value=0.1439)
params.add('V', value=2.561)
params.add('Tlag', value=0.8454)
params.add('E0', value=350)
params.add('a', value=0.005)
params.add('b', value=0.15)
params.add('g', value=0.001)
params.add('AaDose', value=1)  # Example dose, adjust as needed

# Time points for simulation
t_span = (0, 24)
t_eval = np.linspace(0, 24, 100)

# Initial conditions
y0 = [params['AaDose'].value, 0, params['E0'].value]

# Solve the differential equations for initial conditions
sol = solve_ivp(pk_pd_model, t_span, y0, args=(
    params['Ka'].value, 
    params['Ke'].value, 
    params['V'].value, 
    params['Tlag'].value, 
    params['a'].value, 
    params['b'].value, 
    params['E0'].value, 
    params['g'].value, 
    params['AaDose'].value), 
    t_eval=t_eval
)

# Extract results
time = sol.t
E = sol.y[2]

# Data for fitting (replace with your actual data)
noise_level = 5
EObe = E + np.random.normal(0, noise_level, len(E))

# Corrected objective_function
def objective_function(t, Ka, Ke, V, Tlag, a, b, E0, g, AaDose):
    """Calculates the model prediction for a given set of parameters and time points."""
    # Initial conditions
    y0 = [AaDose, 0, E0]  # Use the current values of E0 and AaDose
    # Solve the differential equations
    sol = solve_ivp(pk_pd_model, t_span, y0, args=(Ka, Ke, V, Tlag, a, b, E0, g, AaDose), t_eval=t)
    # Debug: Print the shape of the output
    print(f"Shape of predictions: {sol.y[2].shape}")
    # Return the model prediction for E
    return sol.y[2]  # This should be of shape (len(t),)

# Create model for fitting
model = Model(objective_function, independent_vars=['t'])  # Specify independent variable
result = model.fit(EObe, t=time, params=params)

# Print fitting report
print(result.fit_report())

# Plot results
plt.plot(time, E, label='Simulated E')
plt.plot(time, EObe, label='Observed E (with noise)')
plt.xlabel('Time')
plt.ylabel('Tumor Volume (E)')
plt.legend()
plt.show()

'''

# 2024/10/4 

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Parameters for the PK model  PK模型
Ka = 4.0488   # Absorption rate constant
V = 2.561     # Volume of distribution
Ke = 0.1439   # Elimination rate constant
Tlag = 0.8454 # Lag time

# Parameters for the Tumor Growth Inhibition Model   肿瘤增长/抑制模型与药物效应模型(PD模型)
E0 = 350      # Initial Tumor Volume
a = 0.005     # Growth rate constant
b = 0.15      # Deceleration rate constant
g = 0.001     # Drug decay rate constant

# Initial conditions   肿瘤初始体积
Aa0 = 0       # Initial amount in absorption compartment
A10 = 0       # Initial amount in central compartment
E0 = 350      # Initial tumor volume

# Time points for the simulation
t_eval = np.linspace(0, 24, 100)  # 24 hours divided into 100 time points

def model(t, y, Ka, V, Ke, a, b, g):
    Aa, A1, E = y
    C = A1 / V
    dAa_dt = -Ka * Aa
    dA1_dt = Ka * Aa - Ke * A1
    dE_dt = E * (a - b * np.log(E / E0)) - g * C * E
    return [dAa_dt, dA1_dt, dE_dt]

initial_conditions = [Aa0, A10, E0]

# Solve the differential equations
solution = solve_ivp(
    model, 
    [0, 24], 
    initial_conditions, 
    args=(Ka, V, Ke, a, b, g), 
    t_eval=t_eval
)

# Extracting results
Aa, A1, E = solution.y

# Observed tumor volume with error
EEps = 1  # Error term
EObe = E + EEps

# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(solution.t, EObe, label='Observed Tumor Volume (EObe)', color='b')
plt.xlabel('Time (hours)')
plt.ylabel('Tumor Volume')
plt.title('Tumor Growth Inhibition Model')
plt.legend()
plt.grid(True)
plt.show()


